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Blocking probability is the most frequently required performance charac- 
teristic in traffic studies of complex central office switching networks. Deter- 
mining this quantity without actual measurement is a difficult task. To aid 
the com mimical ions system designer, a simulation program has beenprepared 
lohich produces useful estimates of blocking probability for a large class of 
networks. The program is based on a simplified mathematical model for the 
analysis of switching networks developed by C. Y. Lee, and thus differs 
from conventional simulators in that it simulates a mathematical model 
rather than a traffic-handling system. 

Although Lee's model is widely used, its utility has been limited by com- 
putational difficulties encountered in networks of realistic size and com- 
plexity. This limitation is in most practical cases removed by the program, 
which features rapid input preparation, short computer runs, and specifica- 
tion of the desired precision of the results as an input parameter. Moreover, 
the program allows for the incorporation of more a priori information 
about the actual behaviour of switching networks than is included in Lee's 
model, thereby leading to a more accurate estimate of blocking probability. 

The simulator has been programmed for the IBM 7090 computer, but 
the concepts arc machine independent. 
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Determining the traffic performance of complex multistage central 
office switching systems without actual measurement can be a major 
problem for the communications engineer. While probability theory has 
been successfully applied to a wide variety of telephone traffic prob- 
lems, 123 a precise formulation of a mathematical model completely 
describing the multistage switching system has thus far not been found. 

No systematic approach exists which completely accounts for the 
gross complexity encountered in large-scale congestion systems, but 
several authors have contributed significantly to the theory, notably 
C. Jacobaeus, 46 K. Lundkvist, 5 A. Jensen, 7 C. Y. Lee, 8 A. Elldin, 9 R. 
Fortet, 10 and P. LeGall. 11 More recently, V. E. BeneS 12 " 18 has initiated 
"an attempt to describe a comprehensive point of view towards the 
subject of connecting systems." 12 Although the engineer does not yet 
have a comprehensive theory, he does have a valuable tool in computer 
simulation. 

Simulation of telephone traffic flow has a long history in the Bell 
System. As early as 1907, a rudimentary simulation was undertaken to 
improve switchboard performance. Artificial traffic was generated by a 
card-drawing technique, and the simulation was used to verify a semi- 
mathematical analysis of the loads which could be handled by a team 
of operators meeting an average delay criterion. In the ensuing years, 
simulation techniques have been aids in the study of complex traffic 
problems, such as the effect of limited sources on graded multiple capac- 
ities, the efficiency of random slipped multiples, the capacities of various 
alternate routing plans, and the distribution of delays under various 
trunking plans. The traffic load capacity of the No. 1 crossbar network 
was largely determined by the load-loss relationships in the link and 
junctor patterns obtained from elaborate simulations begun in 1930. 
This was the first time that the capacity of a largely complete system 
had come under study by simulation methods. A 10,000-line No. 5 
crossbar office was simulated in 1948 by a specially designed machine 
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which coordinated the efforts of four operators, providing significant 
data for the traffic engineering of this system.' 6 

In recent years the high-speed electronic digital computer has proved 
an effective tool in large-scale traffic simulations. 17-22 For this class of 
simulations, special computer programs are written which usually con- 
tain : 

(i) a logical description of the system under consideration, 
(ii) a procedure for generating and offering traffic to the system, and 

(Hi) a method for extracting and recording the desired system charac- 
teristics. 

These programs may be called "special-purpose" simulators — in the 
sense that they are written for the purpose of studying a specific traffic- 
handling system. A variety of performance data may be obtained, in- 
cluding: 

(i) probability of blocking at various loads (load-loss data); 
(ii) delay distribution, including average delay on calls delayed; and 

(Hi) mean queue lengths. 
Of these, the most frequently required datum is the probability of loss 
(or delay) . 

These simulation programs have produced a large amount of useful 
information, but their application has not been widespread because of 
the considerable programming effort required. To reduce programming 
effort, various "general-purpose" traffic simulation programs have been 
written.'- 1 2r, ' 2fi However, it must be understood that each program is 
only "general" with respect to a particular class of traffic systems. 

The multistage central office switching system is an example of a class 
of traffic-handling systems for which no general-purpose simulator has 
heretofore been written, although much has been accomplished by 
special-purpose simulations written for specific switching network 
arrangements. 27 Because the use of these network simulation programs 
has been greatly restricted by the cumbersome programming and input 
preparation required, a strong need has developed for a quick, easy-to- 
use, general -purpose simulation technique. 

To meet this need, the authors have developed a computer program 
which, with a minimum of user effort, will produce useful estimates of 
blocking probability for a very large class of multistage switching net- 
works. The program is based on a simplified mathematical model of 
switching networks developed by C\ V. Lee, 8 and differs in approach 
from programs referred to earlier in that a complete description of the 
traffic-handling system itself is not given to the computer; rather, the 
program simulates the behavior of Lee's analytical model. Deriving its 
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name from this view of its operation, the program has been acrony- 
mously named NEASIM — NEtwork Analytical SIMulator. 

While Lee's model is probably the most widely used analytical model 
for multistage matching networks, its utility has been severely limited 
by the computational difficulties associated with networks of realistic 
size and complexity. This limitation is in most practical cases removed 
by the program, which features rapid input preparation and short (hence 
economical) computer runs. Furthermore, unlike many simulations in 
which the reliability of results must be assessed on an a posteriori basis, 
the analytical simulator admits of an a priori appraisal so that desired 
precision becomes an input parameter. 

The probability linear-graph model, basic to the NEASIM approach, 
is described in Section II, where its use in switching network analysis 
is explained. The philosophy of the program together with a general 
description is presented in Section III. Two key program routines are 
described in Section IV. Reliability considerations are given in Section 
V. Finally, Section VI discusses program modifications which can in- 
crease the validity of the NEASIM estimate. 

II. THE PROBABILITY LINEAR-GRAPH MODEL 

In 1955 C. Y. Lee, 8 extending the earlier work of Kittredge and 
Molina, presented a simplified mathematical model for the analysis of 
switching networks. Since the NEASIM program simulates the be- 
havior of this model, a description of the model is given (Section 2.1). 
An example to illustrate how the model is applied is given in Section 
2.2; the computational difficulties which may be encountered in the 
analysis of practical networks are then discussed — thus pointing up the 
need for the simulation program. Section 2.3 explores further the notion 
of blocking probability for a network and the use of the linear-graph 
model in determining this quantity. 

2.1 The Model 

Consider a crosspoint network in which each input can be connected 
to any output by the operation of appropriate crosspoints. Let P t (j,k) 
be the probability that all paths through the network between input j 
and output k are busy at time t. 

Associate with each link U of the network a binary- valued random 
variable X t (i) whose value represents the state of the link at time t. 
The NEASIM convention is: represents idle and 1 represents busy. 



NETWORK ANALYSIS SIMULATION PROGRAM 9G9 

Since the concern of telephone traffic engineers is with "busy-hour" 
traffic, it is usually assumed that 

(a) the busy-idle distributions of the link random variables are sta- 
tionary (or homogeneous) in time. 

That is, for any N, 

Pr{X ln u ' = 8„ ; n = 1, ■•■, V) = Pv{X tn+h (i) = 8 n ; n = 1, •■•, .V) 

for all h, where 8„ = or 1 and the index i runs over all the links in the 
network. A consequence of this assumption is that P t (j,k) is also sta- 
tionary in time, for all j and /,-, and the subscript / will henceforth Ik; 
omitted. 

Let us now fix our attention on two terminals, one on each side of 
a crosspoint network in which 

(b) all of the switches are nonblockiny* 
and in which 

(c) there is no connection path between any switch and itself. 

Then the configuration of possible paths through the network between 
the two terminals can be represented by a two-terminal cycle-free linear 
graph with directed branches, in which the nodes of the graph represent 
network switches and the directed branches of the graph represent 
network links. Consider, for example, the network depicted in Fig. 1 . 
The possible path configuration seen between terminals A and B is 
indicated by heavy lines, and the corresponding graph is as shown in 
Fig. 2(a). 

Next, assume that 

(d) P(j,k) is independent of j and k, 

and we can speak of P(j,k) = B as the probability of blocking of the 
network. The notion of blocking probability will be further explored in 
Section 2.3. 

Finally, Lee makes the simplifying assumption: 

(e) the link random variables, X " , are independent. 

This assumption, which is frequently made to render analysis managea- 
ble, is the principal weakness of the model and will cause the results to 
depart from reality in varying degrees — depending on the particular 
network. In general the model will tend to overestimate blocking. The 
problem of obtaining realistic results is discussed further in Section VI. 
Each of a large class of switching networks can therefore be repre- 



* Lee's requirement that the .switches he nonblocking is actually not restric- 
tive, and can he relaxed by an appropriate adjustment of the link occupancies. 
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SWITCHING NETWORK 

Fig. 1 — Simple crosspoint network with possible paths between terminals A 
and S indicated by heavy lines. 
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Fig. 2— (a) GRAPH of the network of Fig. 1, with occupancy p on each of 
the branches, (b) The subgraphs of the GRAPH. 
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sented by a simplified model called, by Lee, a two-terminal probability 
linear-graph in which 

(i) switches are represented by nodes, and links by directed branches, 
and 

(ii) assumptions (a)-(e) hold. 
The mathematical object, the two-terminal probability linear-graph, 
will be referred to as GRAPH in the sequel. 

Once the GRAPH of a network, is obtained, the calculation of the 
blocking probability can proceed in a straightforward manner. 

2.2 An Example 

As a first illustration, consider the GRAPH of Fig. 2(a). Let E be the 
event that there is a path through the GRAPH,* Ei be the event that 
branch 6, is idle, and p, = Pr{6j is busy) . Then 

E = Ai U .1, U.4 :) U.l, 

where the paths, A , , are 

a , = e x n e, n e, 

,i., = e x n e a n e 8 

a, = e, n e 6 n e 1 

a, = /•;., n A' fi n /•;, . 

The blocking probability is then 
B = \ - ?t{E\ 

= 1 - PrMx U ,1, U,1 ;( U .1,1 

= l - EPrMiJ + EPviAiOAj) 

- £ Pr|,l, f\Ajf) A,\ + PrJ.li D A, C\ A, D A<\, 
ij,k-i 

• <j<k 

Now the assumption of independence gives 

Pr{7? ( ••• Ej\ = PrjA\| ••• Pr{#,} 

= <H • ■ ■ <l, <H = 1 - Pi 



* Having adopted the (1RAPH model, we speak of a "path through the 
CJRAPH" rather than "a path through the network," and say that "a branch is 
busy or idle" rather than "a link is husy or idle." 
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whence, for the case in which p, = p for all i, 

B = q s - 4q 7 + 2r/ + 4g 6 - 4g 3 + 1. (1) 

In general, given a GRAPH G with m different link occupancies 
pi , • ■ • , p m , the procedure just illustrated will yield the blocking 
polynomial of the GRAPH : 

B (! = Baipi , ■ ■ ■ ,Pm). 

As a computational tool the utility of the GRAPH model decreases 
with increasing complexity of the GRAPH'S geometrical structure. 
When the GRAPH geometry grows more complex, the blocking poly- 
nomial B a becomes cumbersome — admitting a greater possibility 
for error in its determination. Moreover, once B Q is found, one still has 
to substitute numerical values for p x , ■ ■ • , p m into the polynomial 
to obtain a result. As an example, for the GRAPH of moderate com- 
plexity shown in Fig. 3, the blocking polynomial is 




OCCUPANCY = p Pz 

Fig. 3 — A GRAPH of moderate complexity. 



where : 
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B = p> 

+ //r/(8,l 4 ) 

+ p 6 q 2 (4A* + •24A i B 2 ) 

+ pq 3 (24A 4 B i + 32A 3 B 3 C) 

+ ? A/(8,lV' 4 + 48,1-BV 2 + SB'D + (>£ 8 ) (2) 

+ p 3 q\Z2AB*C'D + 24fl 4 C 4 ) 

+ pV(4C 8 + 24f?('*D 2 ) 

+ pq 7 (SC i D i ) 

+ r/(7) 8 ) 

■ t = p + gpz 

£ = p + qpS 

C = p + 9p 2 3 

1) = p + qp-i 

3=1- p. 

Faced with computing 5 in ( 2 ) over a range of occupancies, an engineer 
would surely resort to a computer. The essential computational diffi- 
culty is that one is confronted with expressions of the form 

Prj.l, U ••• [)Aj] 

where the paths ^4, , • • • , Aj are not disjoint. 

It is interesting to note that the method of approach just described — 
which may be called the "path enumeration approach" — is not the 
only way to proceed and is indeed not the most efficient. A second pro- 
cedure for finding the blocking polynomial may be called the "combi- 
natorial approach" and is best illustrated by example. 

Since all the branches in the ( i RAPH of Fig. 2(a) arc busy with proba- 
bility p, a moment's reflection shows that the blocking probability can 
be written as 

B = p 1 + qp-B sxlW T AV \> 1 + pq- Subgraph 2 + 7 J -^s„l,„r: 1 ph3 



074 THE BELL SYSTEM TECHNICAL JOURNAL, MAY 1904 

where 5 a ub B « p h i , Subgraph 2 , £sub«r»ph 3 are, respectively, the blocking 
probabilities of the GRAPHs indicated in Fig. 2(b). But by inspection 
we have 

BsubBTQphi = (1 - <h~ 

BsubBraph 2 = (1 — ( l) 
BsubgrnphS = [1 - (1 - P 2 )qf, 

so that 

B - p 2 + 2pq(\ - q-f + '/[I " (1 - P 2 )qf 

is the blocking polynomial ( 1 ) expressed in another form. 

The GRAPH of Fig. 4 is more representative of the type of GRAPH 



Fig. 4 — Typical GRAPH geometry of realistic central office network. 
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encountered in modern centra] office networks. Determination of the 
blocking polynomial in this case — by either the path enumeration or 
combinatorial approaches — is a formidable task indeed. If the GRAPH 
model is to be useful to the modern engineer, some means of handling 
such complex GRAPHs must be made available. 

When seeking performance measures of a network, the engineer 
is not really concerned with the explicit polynomial representation 
BtiiPi , " " • > P>n); what he would like is a curve (or set of curves) dis- 
playing this functional relationship. The simulation program described 
in the following sections, when given the link occupancies pi , ■ • • , p„, , 
produces, with predictable precision, a numerical estimate of B G . 

2.3 Blocking Probability and lire Linear-Graph Model 

When several kinds of traffic arc handled with different disciplines 
in a network which may have several characteristic graphs, blocking 
probability for the network can only be meaningfully defined relative 
to the persons or terminals encountering blocking. For example, the 
blocking encountered by a call originating and terminating in the same 
central office may not be the same as the blocking encountered by an 
incoming call. In general, a network is required to handle several 
"classes" of connections and the engineer is concerned with the blocking 
probability for each of these classes. We shall limit our discussion of 
blocking probability to one class, that is, a subset of all input-output 
pairs in which each input-output pair has the same graph and for which 
it is reasonable to suppose that the traffic between every input-output 
pair is identical* with that of every other pair. Without loss of generality 
we can therefore assume that the network has one class, so that when 
speaking of the blocking probability of the network, we shall mean the 
blocking probability of the class. (These remarks form the basis of 
assumption (d) of Section 2.1. ) 

Having thus limited ourselves to one class of connections, we have still 
to define the blocking probability of a network in a manner which is 
in agreement with the generally familiar definitions. Here we again en- 
counter difficulty. The authors agree wholeheartedly with Benes (Ref. 
15, p. 2805), that "In fact, not even the definition (let alone the calcu- 
lation) of the probability of blocking has received adequate treatment 

Syski (Ref. 3, p. 198), after a long series of prefatory remarks, defines 
two quantities, time congestion S(t) and call congestion ir(l) for the case 



* That is, every input calls every output at the same rate with the same hold- 
ing time distribution and with the same lost calls disposition. 
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of a simple full-access trunk group. The time congestion is the probabil- 
ity that all trunks in the group are busy at time /; the call congestion 
is the conditional probability that the group is blocked when a call 
arrives at the instant t. Under the input assumptions and for equilibrium 
conditions, time and call congestions arc independent of time and are 
denoted by S and t, respectively. 8 is equivalent to the fraction of time 
during which congestion is encountered, while ir is equivalent to the 
fraction of calls encountering congestion, and the two quantities will, in 
general, differ. 

Of these, of course, the measure of most concern to the network de- 
signer is call congestion. Now if it is assumed that calls originate com- 
pletely independent of the state of the network, time congestion will 
equal call congestion. Such an assumption is unjustified if calls cannot 
originate from busy lines, since time congestion conventionally includes 
busy line periods while call congestion excludes them. It is, however, 
reasonable to assume that idle pairs of terminals originate calls at a 
constant rate independent of the state of the network. In particular, 
there must be no change in the calling rate after a blocked call. If, under 
this assumption, time congestion is modified to include only periods in 
which both lines are idle, it will be equal to call congestion. Actually, 
even if the foregoing assumption is not met, the time between calls is 
likely to be much longer than the time taken by the network to return 
to equilibrium, so that, again, the modified time congestion will be close 
to the call congestion. 

With a suitable choice of branch occupancies, Lee's model allows the 
computation of call congestion. Alternatively, the branch occupancies 
may be chosen so as not to reflect the requirement that only idle ter- 
minals are to be considered, thus allowing the computation of time 
congestion. In either case the computed results will be subject to the 
error introduced by inaccurate assumptions in the model. 

Before viewing the probability linear-graph model in the light of the 
above remarks, it is well to make an observation on the underlying 
philosophy of the model. Let us perform the following conceptual experi- 
ment in a real network under a particular set of (equilibrium) traffic 
conditions. Suppose that we fix our attention on a particular representa- 
tive input-output pair (j,k) and examine closely that portion of the 
network seen between them, i.e., their graph. It is reasonable to believe 
that, were we to examine the detailed traffic pattern within the graph 
for a sufficiently long time, Ave would ultimately come to have complete 
knowledge of the busy-idle state behavior of the graph links under the 
particular traffic conditions. The experiment could be repeated under 
other equilibrium traffic conditions, so that we would eventually be able 
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to describe completely the behavior of the graph links under all equilib- 
rium conditions. Assuming that the particular input-output pair and 
connection graph studied are truly representative of the entire network 
(or at least of an entire connection class), we could then "discard" the 
rest of the network and determine the blocking probability of the net- 
work under various equilibrium conditions by computations or simula- 
tions based only on the connection graph and our complete knowledge 
of its behavior. 

That such complete knowledge could be obtained is a practical im- 
possibility. In the absence of complete knowledge, assumptions can be 
and, indeed, must be made about the detailed behavior of the graph 
based on such a priori knowledge as we have. It is reasonable to suppose 
that, as the assumptions made approach the real behavior of the graph, 
the blocking probability determined from the graph will approach the 
real blocking of the network. A belief in the fundamental soundness of 
this reasoning constitutes the basic philosophy of the graph model 
approach to the determination of blocking probability, beginning with 
Molina and continuing with 0. Y. Lee and the NEASIM program. 

The assumption made by Lee of link independence [(e) of Section 
2.1], while obviously omitting much a priori knowledge of graph be- 
havior, possesses the practical advantage of allowing computation of the 
blocking polynomial where graph geometry docs not prohibit. The basic 
NEASIM program allows evaluation of the polynomial for all graphs 
meeting Lee's restrictions. 

The utility of the results obtainable from Lee's model is well known. 
When the specific branch occupancies are chosen rationally,* the calcu- 
lated blocking agrees well enough with real blocking figures (obtained 
from full-scale simulation or measurement) for many engineering and 
design purposes. Accuracy can be improved by "calibrating" Lee's 
results against real values where they are available. Furthermore, com- 
puted values lacking in absolute accuracy will reveal relative differences 
between networks and between various traffic conditions in the same 
network. 

The NEASIM program, to which the remainder of this paper is 
devoted, is basically designed to estimate the value of the blocking 
polynomial. This portion of its design and use is described in Sections 
III-V. The design also allows additional assumptions regarding the 
detailed traffic behavior of connection graph link states to be incor- 
porated in the graph model. When more a priori information is included, 



* Thiil is, cIkpscd to reflect the requirement that the input-output terminals 
j,k are idle by (usually) subtracting the load contributed by the terminals j,k 
from the assumed carried link loads. 
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the validity of the simulation results improves as anticipated. This 
aspect of the NEASIM program is described in Section VI. 

III. THE NEASIM PROGRAM 

This part of the paper is devoted to a presentation of the program. 
The basic viewpoint or philosophy of the NEASIM approach is given 
in Section 3.1. A general description of the program is contained in 
Section 3.2. 

3.1 Philosophy of (he Program 

We saw in the preceding section how the GRAPH model provides — 
in theory at least — a method for the analysis of a large class of switching 
networks and how the model is impractical as a tool for networks with 
complex geometrical structure. To evolve a practical tool, we shall 
change our method of approach. 

Suppose we are given a GRAPH G with branch occupancies 

Pi , ■ ■ ■ , V-n ■ 

Until now we were concerned with the explicit blocking polynomial 
B„ = B a (pi , • • • , p,„). However, we can think of B G (pi , • • • , p m ) as a 
curve in Euclidean (m + 1) -space, and set as our goal a precise ap- 
proximation to this curve. This point of view immediately suggests the 
use of simulation techniques. 

For example, if all the branches in the GRAPH of Fig. 2(a) are busy 
with probability p, then B G = B G (p) is a curve in 2-space. If a com- 
puter simulation were to be performed to give an approximation to 
B a (p), then we would require a computer program containing an 
algorithm which, when repeated n times, will produce an estimate 
B G "\p) such that 

lim B M (p) = B (l (p) 

and for which wc have confidence limits on the absolute error, 

\B G in \p) -B a (p)\, 

for all n. 

Consider the eight branches of the GRAPH of Fig. 2(a). Although 
all of them are busy with probability p, at any one instant each of the 
branches is either busy or idle — and for this particular configuration 
of busy and idle branches there is or is not a path through the graph. 
Thus, in each repetition of the above mentioned algorithm, 
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(*) we assign busy-idle states to each of the eight branches in such 
a way that probability of busy in each case is p; 

(ii) after the assignment has been made we determine whether or 
not there is a path through the graph for the given assignment. 
Let B a M (p) represent the proportion of n repetitions of the algorithm 
when no path through the graph was found. We would then expect that 
lim Bo n) (p) = B a (p); that is, we would expect our estimator to 

n-»oo 

converge to the true blocking curve. In general, if the GRAPH G has 
m different occupancies Pi , ■ ■ ■ , p,„ , the program should produce an 
estimator B a { "\pi , • • • , p m ) converging to B a (p\ , ■ ■ • , p m ). 

The preceding heuristic remarks were intended to outline the essential 
approach taken by the NEASIM program. The remainder of Section 
III is devoted to a description of the program itself. Sections IV and V 
contain the arguments which show that the estimator indeed converges 
to the true blocking curve, and how confidence statements about the 
precision of the results are obtained. 

3.2 Program Description 

The following four sections describe the NEASIM program. Section 
3.2.1 takes a "macroscopic" point of view, beginning with an account of 
the input and then proceeding to give a broad outline of the program. 
Section 3.2.2 sketches the layout of data in the computer memory. 
Section 3.2.3 discusses the organization and operation of the NEASIM 
algorithm. The salient features of the program flow are shown in Fig. 
5. Storage requirements and execution speed are given in Section 3.2.4. 

3.2.1 Macroscopic Description 

NEASIM was written for the IBM 7090 computer. The input con- 
sists of punched cards which we categorize as Graph Definition Cards 
and Simulation Definition Cards. Since the notion of a probability- 
linear graph implies a geometrical configuration together with an oc- 
cupancy assignment on the branches, the computer must be supplied 
with both these types of information. The geometrical configuration is 
read into the computer via the Graph Definition Cards and the various 
occupancies arc read in via the Simulation Definition Cards. 

Graph Definition Cards 

The information punched on these cards includes 
( 1 ) the total number of nodes, 
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Fig. 5 — Flow diagram for the NEASIM program. 
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( 2 ) the total number of branches, 

(3) the total number of branch occupancies, 

(4) the occupancy associated with each of the branches, and 

(5) the interconnection scheme between the various nodes. 

Simulation Definition Cards 

Suppose there are m different branch occupancies for the GRAPH 
in question. These in values are punched on a Simulation Definition 
Card. For every set of values of the occupancies p\ , ■ • ■ , p m that 
is, for every point on the blocking curve B a (pi , •• • , p,„) — there is 
one Simulation Definition Card. 

The estimator B,; U "(p l , ■ ■ • , »,„.) will converge to B„(p\ , ■ • ■ , p,„) 
as n —* co . But the computer must be instructed as to when to terminate 
the run. Now, the NEASIM process is such that it is possible to request 
that the error | /? < " ) — B a | lie within a given percentage of the true 
value B G . This "desired precision" information is also punched on the 
Simulation Definition Card. It may happen, however (as will be the 
case whenever B G is very small) that, in order to obtain the requested 
percentage error, the total number of repetitions of the NEASIM 
algorithm will perforce be exceedingly large. Since a lengthy computer 
run is economically undesirable, and since a high degree of precision is 
usually not needed when the blocking probability is so very low, an 
upper limit on the number of repetitions, n, is also supplied to the com- 
puter by being punched on the Simulation Definition Card. If, after a 
run has been made, a greater degree of precision is still needed, it is 
possible to "pick up where we left off" and continue the simulation in 
another computer run. 

The Program 

The Graph Definition Cards are read into the computer first. With 
this information the program constructs, in effect, a map in the com- 
puter memory of the geometrical configuration of the GRAPH. More- 
over, the program associates with each branch an occupancy p, — whose 
value is as yet unspecified. The first Simulation Definition Card is then 
read in, and the program now assigns the appropriate values to 
Pi , • • ■ , p,„ . Once this information is obtained, the program is ready 
to execute the NEASIM algorithm, which consists essentially of two 
parts : 

(i) a busy-idle assignment is made on all of the branches in accord- 
ance with the specified occupancies p v , • • • , p m ; 
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(n) the presence or absence of a path is determined for the particular 

assignment. 

These two steps are repeated again and again until such time as the 

estimate Ba H) (pi , ••• ,Pm) has been found. 

The results for this point on the blocking curve are printed out, and 
the second Simulation Definition Card (if there is one) is read. The 
program then goes through the preceding steps for this second set of 
values of pi , • • • , p m until the estimate for this point on the blocking 
curve has been obtained. When all the Simulation Definition Cards 
have been processed, the program run ends. 

3.2.2 Memory Organization 

After the Graph Definition Cards have been read in, the program 
prepares five main tables as follows: 

(1) Node table 

(2) Node-Link table 

(3) Branch tables 

(4) Occupancy tables 

(5) Linkage table. 

The Node Table 

The size (i.e., the number of words) of the Node table is equal to the 
number of nodes in the GRAPH, there being a one-to-one correspond- 
ence between the words in this table and the nodes in the GRAPH. 
These words are used by the program to indicate, after a particular 
iteration of the NEASIM algorithm, whether or not there is a path 
from each particular node to the first node. 

The Node-Link Table 

The size of the Node-Link table is also equal to the total number of 
nodes in the GRAPH, with each word corresponding to a particular 
GRAPH node. For each node the table indicates 

(i) the number of branches leaving the node — connecting to nodes 
more distant from the origin, and 

(ii) a reference to a section of the Linkage table where further infor- 
mation on each branch is stored. 

The Node-Link and Linkage tables together constitute the program's 
map of the GRAPH geometry. The other tables provide storage for 
busy-idle indications and path information. 
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The Branch Tables 

There are as many Branch tables as there are different occupancies 
in the GRAPH. For occupancies pi , ■ ■ ■ , p,» there will be m Branch 
tables, which we denote by BRT (1) , • • • , BRT ( „ () . The size of BRT (l) 
equals the number of branches in the GRAPH which are busy with 
probability p< . There is a one-to-one correspondence between a particu- 
lar word in BRT (l) and a particular branch in the GRAPH. The associ- 
ation between the words in the Branch tables and the particular GRAPH 
branches is part of the information stored in the Linkage table. The 
Branch table words are used by the program to store the busy-idle state 
of every GRAPH branch on each iteration of the NEASIM algorithm. 
The total storage required for the Branch tables equals the number of 
branches in the GRAPH. 

The Occupancy Tables 

There are in Occupancy tables, OCT ( d , • • • , OCT (m) , corresponding 
respectively to BRT(i) , • • • , BRT ( ,„) . The size of OCT( t) is an input 
parameter of the program chosen to be large compared with BRT (() . 
Every bit in the table OCT (l) contains a binary one with probability 
Pi . The Occupancy tables are used by the program to supply random 
busy-idle states for assignment to the branches of the GRAPH on each 
iteration of the NEASIM algorithm. 

The Linkage Table 

As previously mentioned, this table contains the detailed interconnec- 
tion information between the various nodes in the GRAPH. The size 
of the table is equal to the number of branches in the GRAPH. Each 
word in the table represents a branch, say bj , in the GRAPH and 
contains 

(i) the address of a word in the Node table which corresponds to 
the node to which bj leads, and 

(ii) the address of a word in a Branch table which stores the current 
busy-idle state of the branch bj . 

Consider, for example, the GRAPH of Fig. 2(a) and suppose that 
the occupancy of branches hi and b- is pi ; the occupancy of branches 
bj , bi , b b and /> 6 is p 2 ; and the occupancy of branches 67 and bs is p 3 . 
The tables which the program would prepare are indicated in Fig. 6. 
The symbolic addresses, such as NW1, BRW3, etc., have been chosen 
for illustrative purposes only. 
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Fig. (i — Computer memory layout for the GRAPH of Fig. 2 with occupancies 
pi , p> and p* . 

3.2.3 Organization of the NEASTM Algorithm 

In a previous section it was mentioned that the NEASIM algorithm 
requires 

(i) an assignment of busy-idle states on the branches, and 
(it) a method of searching for the existence of a path. 
These tasks are performed by three program segments which we shall 
call 

(i) the PROBABILITY GENERATOR, 
(it) the BUSY-IDLE ASSIGNMENT, and 
(Hi) the MATCH routine. 
For the sake of program efficiency (speed) the task of assigning busy- 
idle states to the branches is divided into two parts. The function of 
the PROBABILITY GENERATOR is to generate tables of busy-idle 
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hits the Occupancy tables. The BUSY-IDLE ASSIGNMENT rou- 
tine assigns busy-idle states to branches on each iteration of the algo- 
rithm. The function of the MATCH routine is to find whether or not 
there is a path, given a particular assignment of busy-idle states on the 
branches. Discussion of the internal logic of these programs will be 
deferred to Section IV. For the present we will be concerned only with 
their functions. 

After the Graph Definition Cards have been read, enough information 
is available for the program to reserve appropriate space for the Node, 
Node- Link, Branch, Occupancy and Linkage tables. Once space allot- 
ment has been made and the necessary entries filled into the Node-Link 
and Linkage tables, the program reads the first Simulation Definition 
Card. Among other parameters, this card specifies the m different branch 
occupancies desired. At this point the PROBABILITY GENERATOR 
(ills in the Occupancy tables. When this routine has completed its task, 
each bit in OCT (1 ) , i = 1, • • • , in, will contain a binary one with proba- 
bility p, and a zero with probability 1 — pi . ( Recall that the NEASLNI 
convention is thai one represents busy and zero represents idle. ) The con- 
tents of the Occupancy tables remain unaltered for the entire time that 
the program is seeking an estimate for one point on the blocking curve. 

The storage of many computers is organized into sequentially num- 
bered words, each of which consists of a fixed number of contiguous 
bits, and each of which is addressable by the stored program for logical 
and algebraic manipulation. The number of bits in an IBM 7090 word 
is 30, so that each word in the Occupancy tables represents 30 independ- 
ent busy-idle states — thus allowing for 30 independent repetitions of 
the NEASIM algorithm. 

After the Occupancy tables have been prepared, various counters are 
set to zero and the Node table is initialized. As the contents of the Node 
table indicate for each node the presence or absence of a path from the 
particular node to the origin, the program takes the attitude of the man 
from Missouri and assumes there is no path until one is proven to exist. 
Thus the words of the Node table are initially set to all 1 's except for 
the first node, which always contains zeros. 

At this point the program enters the BUSY-IDLE ASSIGNMENT 
routine. To assign busy-idle states to the branches, this routine steps 
through each of the Branch tables and for each word in BRT(o , a 
word from OCT (l) is selected at random and its contents duplicated in 
the Branch table word. When BUSY-IDLE ASSIGNMENT has been 
completed, the program enters the MATCH routine. 

Using the linkage information stored in the Node-Link and Linkage 
tables, the MATCH routine performs its logic on the Branch and Node 
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tables to find whether or not there is a path through the graph for each 
of the 36 independent busy-idle configurations. The operation of this 
routine is analyzed in Section 4.2. 

When MATCH has completed its work, the number of l's in the word 
of the Node table corresponding to the terminal node in the GRAPH 
(NW6 in Fig. 6) will be the number of times there was no path through 
the graph in the 36 trials. This number of blocked attempts is noted and 
36 is entered into a "run-length" counter. The Node table is reinitialized, 
busy-idle states reassigned and the algorithm repeated again and again. 
The repetitions will terminate when one of two situations occurs at the 
end of the MATCH routine. Either 

(i) enough attempts have been scored to guarantee the precision 
called for on the Simulation Definition Card, or 

(it) the maximum number of attempts specified on the card has been 
exceeded . 

When the repetitions terminate, the proportion of blocked attempts 
is calculated and the results are printed out. The next Simulation Defini- 
tion Card is read and the process starts over for the next point on the 
blocking curve. 

3.2.4 Storage Requirements and Execution Speed 

The largest GRAPH that can be handled by the present version of 
NEASIM is determined by the core storage available in a particular 
computer. Total storage required is given by the expression 

P + 2(.V + B) + 

where: 

P is the number of storage locations required by the NEASIM 
program itself (about 2000 words), 

iV is the number of nodes in the GRAPH, 

B is the number of branches in the GRAPH, and 

is the storage required for Occupancy tables — typically m X 1024 
where 1 ^ m ^ 8 is the number of occupancy tables. 

The NEASIM algorithm is designed for rapid execution on the IBM 
7090. Average speed depends principally on the size of the GRAPH. 
Typical speeds range from about 600 trials per second (550-branch 
GRAPH) to about 5000 trials per second (48-branch GRAPH). 

IV. PROBABILITY GENERATOR AND MATCH ROUTINES 

The functions of the two routines called PROBABILITY GENERA- 
TOR and MATCH were mentioned in the previous section. The present 
section is concerned with the internal logic of these programs. 
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The algorithm used by PROBABILITY GENERATOR is clue to 
W. C. Jones. 28 It is felt that this heretofore unpublished algorithm is of 
sufficiently widespread interest to be included in this paper. 

A word on notation : throughout Section IV we shall use the notation 
C[A\ to mean "the contents of .4," where A represents some bit in the 
computer memory. (The reader is therefore cautioned to distinguish 
between "bit A" and its contents C[.4].) 

Two computer instructions which will be referred to repeatedly are 
the logical (inclusive) "OR" and logical "AND" instructions. Instruct- 
ing the computer to perform an OR on two bits will guarantee that the 
result will be 1 if, and only if, either or both of the bits contain 1 ; while 
an AND yields 1 if, and only if, both bits contain 1. When we OR bit A 
to bit B, then we shall say that we "perform [A] OR [B] ;" when we AND 
bit A to bit B, then we say that we "perform [.4] AND [B[." 

i.i The I 1 ROB ABILITY GENERATOR 

The purpose of PROBABILITY GENERATOR is to generate the 
occupancies Pi , ■ ■ ■ , p,„ • We mentioned in Section 3.2.2 that this 
subroutine will fill up OCT (l) (i = 1, • • • , m) with 36-bit words each 
of whose bits contains 1 with probability p,-(i = 1, • • • ,m). In the 
sequel we will focus our attention on (a representative) one of these 
30 bits ■ — - keeping in mind that the program is actually working on 
30 bits independently and in parallel. 

Suppose we wish to generate a random binary variable which takes 
the value 1 with probability p(0 < p < 1 ) and to place our result in 
bit X. The algorithm, when terminated, will give PrjC'LY] = 1) = p. 

The first action taken by PROBABILITY GENERATOR is to ex- 
press p as a binary fraction to 10 places.* This fraction is then scanned 
from right to left until the first bit is found which contains a 1. The 
algorithm uses this binary fraction of n ^ 10 places determined by the 
scan. Wc can therefore, without loss of generality, express p as 

p = O.&A-i • • • &2&, bi = 1 ; b 2 , ■ ■ ■ , b n = or 1. 

A digit selected from a random binary number will be referred to as a 
"random bit" in the following algorithm. In a random binary number, 
the value of each digit is 1 with probability h. 

Algorithm: 

(i) Set j = \. Generate a random bit, say 1\ , and store it in X. 
Thus C[X] = r, . 

* The number of places is arbitrary. Ten was chosen to make the round-off 
error smaller than (I.OOl, since occupancies are specified on the Simulation Defini- 
tion Cards to three decimal places. 
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(n) If./ = n, go to step (vii). Otherwise, increase./ by 1 and con- 
tinue. 

(Hi) Generate another random bit rj and store it temporarily in some 
bit, say R. Thus C[R] = r y . 

(iv) If bj = 0, go to step (V). If bj = 1, go to step (vi). 
(v) Perform [R] AND [X]; store in X. Go to step (tt). 
(vi) Perform [R] OR [X]; store in X. Go to step (it). 
(vii) Stop. 

We observe that step (ii) is performed (iterated) exactly n times. 
Let Pj be the probability that C[X] = 1 after the jth iteration. The 
assertion is that P„ = p — i.e., after the algorithm is terminated, X 
will contain 1 with probability p. 

I 'roof of Algorithm: 

Consider two bits .4 and B, and let p A = Pr{C[A] =1), and p B = 
Vv[C[B\ = 1). When the contents of A and B are independent, if we 
perform [A] AND [B], the probability of the result being 1 is 

PaPu , 
while if we perform [A] OR [B], the probability of the result being 1 is 

Pa + Pb - PaPb ■ 
Now, ?v[rj = II = Prjr,- = 0\ = h j = 1, • • • , n. 
Therefore if step (») is executed, 

F; + i=£P; j=h ■••,»- 1 (3) 

while if step (w) is executed, 

P i+1 = i + P, -if ,- i - 1, -•■ ,»- 1 ... 

(4) 

— iP -4- i 

But step (v) is executed only if b j+l = 0, and step (vi) is executed only 
if 6j+i = 1. Hence (3) and (4) can be combined into 

Py+i = 1P/+ I6/+1 J - 1, ■■-,■»- 1. 

Since i'l = £, it follows by induction that 

2 2 2 2 n_I 2" 

and our assertion is proved. 
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4.2 The MATCH Routine 

The .MATCH routine is entered by the program after the busy-idle 
states have been assigned to the branches. Its purpose is to find whether 
or not there is a path through the graph for any particular assignment. 
In the present section we first derive a recursive set-theoretic formula 
which may be used to determine whether there is a path. The remainder 
of the section shows how the recursion is carried out by the computer 
to produce an estimate which converges to the blocking probability of 
the GRAPH. 

In Section 2.2 we saw how the "path enumeration approach" could, 
in theory at least, be employed in determining the blocking probability. 
We again take the path enumeration approach, but from a slightly 
altered point of view: 

Instead of trying to find whether there is a path through the entire 
graph at one fell swoop (which amounts to finding whether there is a 
path from the first node to the last node), we shall try to find whether 
there is a path from the first node to each of the nodes in the GRAPH. 
While this approach may appear to inject an unnecessary complica- 
tion, we will see how, by stepping through the GRAPH in an orderly 
fashion, this approach lends itself naturally to computer programming. 
We begin by investigating, somewhat further, the geometrical structure 
of a GRAPH. 

Consider the general GRAPH (whose structure is shown schematically 
in Fig. 7) and suppose that there arc a total of v nodes. Each GRAPH 
has a first node, N\ , a last node, N v , and several intermediate "stages" 
of nodes. The notion of "stage" is made more precise by the following 
definition: a node N is said to be in stage s(s = 1,2, • • • ) if , and only if, 
all paths from Ni to N contain no more than s — 1 branches, and there 
exists at least one path from Ni to A 7 which contains exactly s — 1 
branches. 

Any GRAPH will thus contain some number S ^ 2 of stages, where 
the first and last ( »S'th stage) consist, respectively, of the single nodes 
Ni and A/, . Let n s be the number of nodes in stage s, s = ],■•■, S (thus 
ih = n s =1); and define a quantity m» by 



/// 



. = £ * ■ 



We choose to order the nodes of the GRAPH in the following manner: 
to each of the n, nodes in stage s (s = 2, • • • , S) assign arbitrarily, but 
uniquely, one of the integers 

w,_i + 1, //?„_! + 2, ■ ■ • , m s _i -J- n„ s = 2, — , S. 
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N 2 




S=l 5 = 2 S = 3 S=S-I S = S 

Fig. 7 — A general GRAPH. 

The nodes of the GRAPH are hence totally ordered and may therefore 
be denoted by Ni , where i = 1, ■ ■ ■ , v. 

When Ni and Nj are connected by one or more parallel branches, 
where JV,- , Nj are arbitrary nodes, we say that, for i < j, Ni is directly 
connected to Nj through each of these branches and write Ni — * Nj . 
We shall suppose, without loss of generality, that if Ni — » Nj, then 
there is only one connecting branch, b/. (This restriction is assumed in 
order to avoid introducing yet another index, but will not affect the 
subsequent results. When a computer instruction involving 6/ is de- 
scribed, we shall understand that the computer is to execute this in- 
struction for all the b/ for which Ni -*Nj.) 

Let us now examine the graph, asking the question for every node 
Nj , "Is there an available path from Ni to Nj ?" Our objective is to 
answer the question for j = v. 

For each node Nj(j > 1 ) consider the branches b/ for which JV< —>Nj. 
Clearly, for each of these branches, if there was no available path from 
Ni to Ni or there is no available path through the branch b/, or both, 
then, and only then, will there be no available path from Ni to Nj which 
passes through Ni . More precisely, let 

Xj be the event: there is no available path from N\ to Nj , 
Yi be the event: there is no available path from N\ to N., 

which passes through Ni , 
Bt be the event: branch b { } is not available (busy). 
Then 
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IV = Xi U B/ 

Xj = n y/ = n (Xi u b/). 

i<i i<i 

Consider all the nodes AT,- such that AT,- — » AT,- . Suppose there are /c, 
such nodes. Call them N %1 , N,- 2 , • • • , A r ,- fc . , where i\ < i> < • ■ ■ < 4- . 
Next, define an event Zj recursively by 

Zj (m) = (X im U Bj) D Z/ m_1) m = 2, • • • , fc,- (5) 

Ni m -> AT,- 

J = 2, • • • , r 

and 

Z/ l) = X tl U B,/ JVi, -► JVi (6) 

./ = 2, • • ■ , v. 

It then follows that 

Z/ k,) = Xj j = 2, • • • , v (7) 



and in particular 



z < fc -> = a; 



where X, is the event "there is no path through the graph." The MATCH 
routine is based upon these formulas. 

We now focus our attention on (a representative) one bit in each 
word of the Node and Branch tables — again keeping in mind that the 
program is actually working on a full word of bits independently and in 
parallel. The program will thus execute 36 simultaneous iterations of 
the MATCH algorithm, which we now proceed to evolve. 

When MATCH is entered, the following state of affairs prevails: 

(i) To each node Af, and to each branch b, J there has been uniquely 
assigned one bit of computer memory, so that we can henceforth refer 
to these bits as bit A^ and bit 6, J . 

(ii) C[b/] = 1 with the appropriate occupancy p/ = Pr| branch 
b/ is busy) • (These p/ were called pi , • • • , p m in Sections II and III.) 

(Hi) The linkage information between the node bits Ni and the 
branch bits b/ is stored in the Node-Link and Linkage tables. 

It follows that representation of the event Zj = X tl U B^' of for- 
mula (6) can be achieved in the computer by performing 

[A',J OR [b h j ] 
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provided that we always set C[N,] = 0; and the event Z/ fc/) = Xj of 
(7) can be represented by executing the following steps. 
(0 Perform [N i{ ] OR (&,/]; store in Nj . 
(ii) Set m = 1. 

(Hi) If m = kj go to step (vi). Otherwise, continue. 
(iv) Increase m by 1. 
(v) Perform ([JV.-J OR [bj]) AND [N } ]; store in Nj . Go to step 

(Hi). 
(vi) Stop. 
We observe that the algorithm may be condensed if we initially set 
C[Nj\ = l(j^2): 
(i) Setm = 1. 

(ii) Perform ([N im ] OR [bj]) AND [A/,]; store in Nj . 
(fit) If m = kj , go to step (v). Otherwise continue. 
(iv) Increase ra by 1. Go to step (ii). 
(v) Stop. 
This initialization of the node bits — C[Ni] = 0, C[Nj] = 1 for j = 
2, 3, • • • , v — was mentioned in Section 3.2.3 and can be interpreted as 
meaning "there is always a path from NitoNi] there is no path from 
Ni to Nj (j > 1) until proven otherwise." 

Summarizing, the steps that the computer may take to represent the 
event Z,*'' = X v are 

(Ml) Set C[Nt] = 0. Set C[Nj] = 1, ./ > 1. 

(M2) For each bit Nj , j = 2, • • • , v, and in the order 

N* , N 3 , ■ • • , N, 

find all bits Ni for which Ni —> Nj and perform 

([Ni] OR [&/]) AND [N } ]; store in Nj . 

Now, the statement "find all bits Ni for which AT, -> AT/' implies 
that the input to the program is such that it specifies to the computer 
which nodes Ni are directly connected to Nj . It was felt that a more 
straightforward input format would be to specify to which nodes, Nj , 
each node N { connects. With the linkage information stored in this 
fashion,* step (M2) may be replaced by the equivalent step 

(M2') For each bit Ni , i = 1, ■ ■ ■ , v — 1, and in the order 

Ni , N 2 , ■ ■ ■ , A 7 ",-! 



linkag 



A careful perusal of Section 3.2.2 will show that this is indeed the way the 
age information is stored in the Node-Link and Linkage tables. 
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find all bits .V, for which .V, — » :V> and perform 

([AM OR [bi 'I ) AND [.V,]; store in Nj . 

Steps (Ml) and (M2') constitute the MATCH algorithm for deter- 
mining whether there is a path through the graph. When the MATCH 
algorithm is terminated, C[N V ] = 1 if, and only if, there is no path, so 
that PrjCLV,,] = lj = Pr{no path}. 

Suppose, the algorithm is iterated n times. In any one iteration the 
event ('[!>/] = 1 is generated with probability p/ and independently of 
the event C[bi] — 1 in any other iteration. But Pr{C[2VJ = 1| after a 
given iteration is clearly only a function of the probabilities Pr{C[6/] = 
1(. Hence, A, will contain 1 with the same probability, B, after each 
iteration, and we conclude that n iterations of the algorithm constitute 
a sequence of /; Bernoulli trials with probability B of success on each 
trial. Let £, be a random variable such that 

& = 1 if C[N,] = 1 after the tth trial 

= if C[N,] = after the ith trial 

and let 



B„ = fc + fc + • • • + £„ . 



B„ = & + h + • • • + £« • 
An application of the Law of Large Numbers (Ref. 29, p. 189) gives 

lim (B n /n) = B. 



Identifying B„ n with /?,,' "'(/>,, ••• , p m ) of Section 3.1 and B with 
B,;(pi , • • • , p m ) of the same section, this last result is equivalent to, 
and proves the assertion that, 

lim Ba"\pi , • • ■ , />,„) = Bo(Pl , • • • , Pm). 

II -»x 

Repeating the experiment often enough and dividing the number of 
times C[N V ] = 1 (i.e., the number of blocked attempts) by the number 
of iterations n (i.e., the number of attempts) will therefore yield a 
precise estimate of the blocking probability of the GRAPH. How- 
large n has to be in order to attain any given degree of precision is dis- 
cussed in the next section. 

V. RELIABILITY CONSIDERATIONS 

Since any simulation is nothing but an experimental measurement, 
and hence subject to statistical fluctuations, it is necessary to assess the 
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reliability of the results. Unfortunately, most of the testing for a given 
simulation must he done on an a -posteriori basis and there is, in general, 
rarely sufficient a priori information on which to base the decision of how 
long the run is to be. The need for such information becomes even more 
important when the cost factor in computer simulations is considered. 
It will be shown below that the NEASIM procedure is one which allows 
for a priori determination of run length. The decision of how long the 
run is to be is based on the desired precision and is made by the com- 
puter. 

In Section 4.2 we introduced the random variable B n , the number of 
times C[N y ) = 1 in n repetitions of the experiment. B n will now be 
called "the number of calls blocked in a simulation run of n calls."* 
The estimator B„/n was seen to converge to the blocking probability 
B, and is indeed a maximum likelihood, mean-unbiased, minimum 
variance estimator. 

To generate the random binary numbers of Section 4.1, NEASIM 
uses the well-known multiplicative congruential method. This method 
for pseudo-random number generation has been discussed, tested, and 
used by numerous investigators 31 " 36 since it was first proposed by 
Lehmer. 37 Although the method has been demonstrated to generate 
35-bit random binary numbers, there is some measure of cyclic behavior 
in the low-order bits. NEASIM forms a word of 36 random bits by 
combining the most random halves (18 high-order bits) of two 35-bit 
words generated by this method. A further check on the randomness is 
provided by NEASIM itself, which, as part of its output, prints the 
generated branch occupancies. We therefore assume that the events 
C[N V ] = 1 are independent for individual trials of the experiment. The 
results of any given computer run should thus be binomially distributed 
(see also Section 4.2), and hence asymptotically normal. 

As an illustration, we again turn to the GRAPH of Fig. 3. For a 
branch occupancy of p = p 2 = 0.5, formula (2) yields B = 0.0525. 
The NEASIM program, for a run of n = 201,600 calls, gave B = 
0.0529 — an error within 1 per cent. (This run took some 40 seconds of 

* The new nomenclature is chosen to indicate a different interpretation of what 
NEASIM is doing: NEASIM looks at the configuration of possible paths through 
the network between two subscribers, takes "snapshots" of the current busy-idle 
states of the links in these paths, and then finds if there is a path for each snap- 
shot. Thus the program is essentially running calls through a portion of the 
network — the portion of the network being a representative one, and hence one 
from which significant statistical data can be extracted to describe the perform- 
ance of the entire network. A similar approach was taken by A. Feiner, W. C. 
Jones, and others, 30 who wrote "abbreviated" simulations in which the busy- 
idle state data were taken from previously run full-scale simulations, but for 
which a new program had to be written for each new network to be simulated. 
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computer time.) To study the normality, the number of blocked calls 
was noted for every 1008 calls run. Since B n is the number of blocked 
calls in a run of n calls, the binomial distribution function in this case is 

Pr{£„ ^ /3) = g ( 10 A 08 ) B\l - B) ms - k 

and the approximating normal distribution function (Kof. 29, p. 172) 
is<I>(.i^ + i), where 

x t = (/ - 1008/?)// 

// = |1008#(1 - B)]~ l 



and 



'I' 



1 r 
■\Z2tt J- a 



e dy 



is the standard normal distribution function. Fig. 8 is a plot on proba- 
bility paper of the cumulative frequency distribution determined by the 
computer, and the theoretical normal distribution line for B = 0.0525. 
In order to obtain meaningful results, it was felt that the run-length 
should be determined by the following criterion: the number of trials 
shall be large enough so as to give 95 per cent confidence that the esti- 
mator lies within a fixed percentage of the true value of B. Since the 
blocking probability is unknown at the outset, this criterion is more 
useful than requiring the estimator to lie in a fixed interval about B. 
Thus, we wish to choose n large enough so that 
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Fig. 8 — Cumulative distribution for GRAPH of Fig. 3 (p = p 2 = 0.5) in a 
run of 201,000 chills in which the number of blocked calls wus noted for every 
1,008 calls run. The points are experimental data; the line is the theoretical dis- 
tribution, *(£0 + j) 
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Pr 



^-B 



n 



< aB) > 0.95 



where a is an input parameter to the program representing the desired 
percentage. Rearranging terms, we get 

B n - nB 



It 



VnB(l - B) 



< a 



nB 



> 0.95 



1 - B\ 

where, for n sufficiently large, the distribution of the random variable 

B n - nB 
VnBQ. - B) 

approaches the standard normal. Our requirement will therefore be 

satisfied if 



or 



- \/^b * " 



> 3.84/1 
a 2 \B 



96 



(8) 



It follows that the number of trials should be increased as the blocking 
probability decreases in order to maintain 95 per cent confidence that 
the estimator lies within a fixed percentage of the true value — as 
intuition dictates. On the other hand, for a fixed B, as n increases, the 
95 per cent confidence limits will narrow. This is displayed in Fig. 9, 
where the program results for the GRAPH of Fig. 3, for B = 0.0525, 
are seen to lie within the confidence limits. 

Furthermore, since (1/B) - 1 ^ ( 1/B) for all B in the unit interval, 
the requirement will be met if 



n > 



3.84 1 



a- 



B' 



But for large n, B tt (B n /n). Hence, making the substitution we obtain 

B „s 3 4 4 . 

a 2 

As long as the number of simulated calls is large enough so that the 
number of blocked calls is at least 3.84/a 2 , there is 95 per cent confidence 
that B n /n is within aB of B. 

NEASIM was written to accept several values of a ranging from 0.05 
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Fig. ( J — NEASIM blocking est i mates for GRAPH of Fig. 15 vs run length 
{,, = v ., = 0.5, B = 0.0525). 

to 1.00. From the specified value, it determines the minimum number of 
blocked calls necessary to guarantee this precision. Fig. 10 shows the 
results of the simulation of the GRAPH of Fig. 3 for a = 0.10 and 0.50. 
In GRAPHs where B is very small, it is usually unnecessary to estimate 
B with a very high degree of precision — especially at the expense of 
costly computer runs. An upper limit for n is therefore also specified. 
The computer then proceeds with the simulation until it either exceeds 
the lower limit on B„ or the upper limit on //. 

In the latter case, the reliability can be assessed as follows. Since 

B„ - nB 



IV 



VnB{l - B) 



< 1.96 > > 0.95 



we can rearrange terms and obtain 

Vv\(IB 2 + eB+f £ 0] ^ 0.95 



where 
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d = n - + 3.8416 n 

e = -(2nB n + 3.8416 n) 

f = b,; 1 . 

It is easily verified that for all n, B n > 0, the parabola dB 1 -f- eB + / 
is concave upward, has real roots, and that B„/n lies between the roots 
whenever B„ < n. Thus for any simulation run, the 95 per cent confi- 
dence interval can be obtained by solving for the roots. 
Now in all cases of interest, the product of the roots, 

f/d = Bn/{n + 3.8416 n), 

can be closely approximated by (B n /n) 2 , so that B„/n can be taken as 
the geometric mean of the roots. Suppose the higher root is B x and the 
lower root B> ; then 
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B x = kBJn 

B, = (l/k)B n /n. 

The roots wore calculated over a range of values of B n and n, and the 
results are displayed in Fig. 1 1 , where k is plotted as a function of B n /n 
for various values of n. As an example, suppose that in a run of 80,000 
calls, B n = 80 were found blocked. Then B n /n = 0.001 and Fig. 11 
gives A- = 1.244. Hence B x = 1.244 (0.001) = 0.00124, B 2 = 0.001/ 
1.244 = 0.000803, and there is 95 per cent confidence that 0.000803 ^ 
B ^ 0.001244. The importance of the case under consideration would 
then determine whether a longer simulation is necessary. 

Since the GRAPH geometry may be as complex as indicated in Pig. 4, 
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it is important to be able to verify that no error was made in mapping 
the geometry into the computer memory. To this end, a program was 
prepared by Miss D. Logan which, using the Graph Definition Cards 
as input and an SC-4020 Microfilm Recorder, draws a picture showing 
the GRAPH geometry. 

VI. POSSIBILITIES FOR INCREASED REALISM 

In the preceding section it was shown that the NEASIM program 
produces precise, reliable results in comparison with Lee's probability 
linear-graph model of switching networks. Unfortunately, in complex 
switching networks, substantial differences may exist between the 
estimates obtained with Lee's analytical technique and actual perform- 
ance determined by field measurement or full-scale (complete) simula- 
tion. These differences appear to be largely attributable to the unreality 
of the assumption [(e) of Section 2.1] that the link random variables are 
independent. 

The dependence that exists between the different links in a network 
stage and between the links in adjacent stages is incompletely under- 
stood, but nonetheless real. Attempts to take account of interlink de- 
pendence by judicious modification of link occupancy values have been 
moderately successful in calculations and may, of course, be employed 
in NEASIM runs. However, the nature of the NEASIM process suggests 
alternate approaches which may ultimately prove to be more fruitful. 
While considerable success has been obtained with some of the tech- 
niques discussed in this section, much remains to be accomplished. 

6.1 Dependence Effects within a Switching Stage 

Two possibilities for increased realism within the confines of a single 
link stage appear worthy of mention. The NF]ASIM program typically 
assigns link stage busy-idle states from a binomial distribution. An 
obvious suggestion (but one upon which little work has been done) 
would be to modify the PROBABILITY GENERATOR and /or 
BUSY-IDLE ASSIGNMENT routines in such a way as to produce 
busy-idle state assignments taken from various distributions — Ja- 
cobacus' E distribution, 40 for example. 

A second approach, which has been extensively used, is to incorporate 
program routines between BUSY-IDLE ASSIGNMENT and execution 
of the MATCH routine. These routines examine the random busy-idle 
assignments and make appropriate assignment changes where GRAPH 
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geometry or other considerations indicate. An example of such a routine 
is one which is designed to insure that there exists at least one "sure- 
idle" branch out of an n X n input switch. If calls are only placed be- 
tween idle network terminals, at least one branch out of an n X n input 
switch must be idle. But if the branches leaving the initial GRAPH 
node are specified at occupancy p, then the program would normally 
make all these branches busy with probability p". The SURE-IDLE 
routine, upon finding such an assignment, will select an initial branch at 
random and make it idle. This action, of course, disturbs the otherwise 
binomial distribution of branch states and should be taken into account 
when branch occupancy values are specified. 

G.2 Interstage Dependence Effects 

An interesting technique for introducing interstage dependence exists 
within the NEASIM framework and is based upon an obvious extension 
of the SURE-IDLE routine just described. Briefly stated, a routine 
could be designed to examine the busy-idle assignments made on the 
input and output branches of each node of the GRAPH. Acting on 
knowledge of the switch geometry and other factors, the "dependencing" 
routine could change the initial busy-idle assignments where necessary 
to make them more realistic. Only a very simple and admittedly in- 
accurate routine has been used to date with, however, a remarkable 
increase in the "realism" of the results obtained. 

The simple-minded DEPENDENCING routine currently employed 
assumes that every GRAPH node is in reality an n X n switch. It 
observes that in an n X n switch each input branch has probability l/n 
of being connected to any particular output branch. It attempts to 
implement its idea of reality by, for each input branch, examining each 
output branch and forcing the output branch state to agree with the 
input branch state with probability l/n. While this routine can produce 
rather quaint effects, such as duplicating the state of one input branch 
on all output branches, it does possess the virtues of simplicity (rapid 
execution) and a basically correct notion of interstage dependence. 

That the use of the SURE-IDLE and DEPENDENCING routines 
can be effective is demonstrated in Fig. 12. The results shown in the 
figure are for a moderately complex eight-stage switching network whose 
GRAPH has 232 branches. NEASIM results with and without de- 
pendencing and sure-idles are compared with the results of a full-scale 
simulation. The improvement in realism possible with the rudimentary 
routines just described appears quite dramatic. 



1002 



THE BELL SYSTEM TECHNICAL JOURNAL, MAY 1964 



0.4 


• 


NEASIM DATA FOR BINOMIAL 








0.2 

> 



A 


NEASIM DATA FOR DEPENDENCING 
AND SURE-IDLE ROUTINES 
FULL-SCALE SIMULATION DATA 




m 












s 


S 


m 0.08 

o 

<r 

■■ 0.06 

U) 

? 0.05 

s: 

O 0.04 

_l 

10 

0.03 
0.02 

0.01 








s 




<T 












^ 












^ .' 


y 










J 


f s 












// 


r 










^ 


/ 










t 


/ 













0.34 



0.36 



0.38 



0.40 0.42 

LOAD 



0.44 



0.46 



0.48 



Fig. 12 — Comparison of NEASIM results (with and without dependencing 
and sure-idles) with full-scale simulation data for a realistically complex network. 



VII. CONCLUSION 

The NEASIM approach to switching network simulation has achieved 
its major goal, the mechanization of a widely employed — if somewhat 
unrealistic — technique for load-loss analysis. The applications of the 
GRAPH model need no longer be restricted by the overpowering com- 
putational difficulties brought on by GRAPH complexity. 

Furthermore, the simulation of an analytical model concept basic to 
NEASIM seems to open new areas for profitable exploration in the 
analysis of switching systems — and perhaps other stochastic systems 
as well. The success of the elementary realism -injecting routines suggests 
that further research along this line may be rewarding. 

Finally, the degree of realism in results attained so far, coupled with 
the ease of application, has produced what amounts to a new tool for 
use in both the design and engineering of new switching networks. It is 
now feasible to achieve relatively complete and accurate load-loss 
engineering data on complex switching systems well in advance of actual 
field experience. 
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